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Using the variational cluster approach (VCA), we study the transition from the antiferromagnetic 
to the superconducting phase of the two-dimensional Hubbard model at zero temperature. Our cal- 
culations are based on a new method to evaluate the VCA grand potential which employs a modified 
Lanczos algorithm and avoids integrations over the real or imaginary frequency axis. Thereby, very 
accurate results are possible for cluster sizes not accessible to full diagonalization. This is important 
for an improved treatment of short-range correlations, including correlations between Cooper pairs 
in particular. We apply this improved method in order to investigate the cluster-size dependence 
of the phase-separation tendency that has been proposed recently on the basis of calculations for 
smaller clusters. While the energy barrier associated with phase separation rapidly decreases with 
increasing cluster size for both hole and electron doping, the extension of the phase-separation re- 
gion behaves differently in the two cases. More specifically, our results suggest that phase separation 
remains persistent in the hole-doped and disappears in the electron-doped case. We also study the 
evolution of the single-particle spectrum as a function of doping and point out the relevance of our 
results for experimental findings in electron and hole-doped materials. 

PACS numbers: 71.10.-w, 74.20.-z, 75.10.-b, 71.30.+h 



I. INTRODUCTION 

Since the discovery of high-temperature superconduc- 
tivity in copper-based transition- metal oxides, a tremen- 
dous effort has been devoted to establish a convincing 
theory that covers the general aspects of their unusual 
and fascinating physics. The attempts are complicated 
by the fact that strong electron correlations play a key 
role in the physics of the cuprates. A central question 
in this context concerns the emergence of small energy 
scales, much smaller than the bare (Coulomb) interac- 
tions between the electrons, which govern the existence 
and the competition of different phases at low temper- 
atures. This can be studied by considering prototypi- 
cal lattice models of strongly correlated electrons. Some 
agreement has been achieved that the relevant physics of 
the cuprate high-temperature superconductors is covered 
by the two-dimensional one-band Hubbard models 

Here Uj denote the hopping matrix elements, n,f is the 
density at site i with spin n% = + n^, \x the 
chemical potential, and U the local Coulomb repulsion. 

In recent years there has been substantial progress in 
the understanding of the ground-state properties of the 
Hubbard model due to the development of quantum- 
cluster theories? 2 - such as cluster extensions of the dy- 
namical mean- field theory (DMFT), i.e. the dynamical 
cluster approximation 3 - and the cellular DMFT— or 
the variational cluster approach (VCA)tf These clus- 
ter calculations confirm the fact that the ground state 



away from half-filling has a non-vanishing supercon- 
ducting order paramete r 8 ' 9 ' 11 ! 12 ) 13 with a pairing in- 
teraction of predominantly d-wave character^ Recent 
VCA calculations&2iiL2 suggest that at low electron and 
hole doping the two-dimensional Hubbard model is in a 
symmetry-broken mixed AF+SC state where both the 
antiferromagnetic (AF) and the superconducting (SC) 
order parameters are finite. This is consistent with re- 
cent cellular DMFT calculations pi 3 - When going to higher 
dopings, the system displays a tendency to phase sepa- 
rate into an AF+SC phase at lower doping and a pure 
SC phase at higher doping. 

The VCA accesses the physics of a lattice model in 
the thermodynamic limit by optimizing trial self-energies 
generated by a reference system. The above-mentioned 
VCA calculations are based on a reference system con- 
sisting of small (2 x 2) isolated clusters tiling the infinite 
lattice. This generates trial self-energies which are very 
short ranged spatially. Hence, there is the obvious ques- 
tion for the robustness of the results as a function of the 
size of the individual clusters. Up to now, it was not 
possible to consider larger cluster sizes and, at the same 
time, reach a sufficient accuracy to resolve the tiny energy 
scale driving phase separation, especially in the electron- 
doped case. The reason is that an accurate evaluation 
of the VCA grand potential has required a full diagonal- 
ization of the cluster Hamiltonian (see Ref. [9| for details) 
which has severely restricted the available cluster sizes. 

The purpose of this paper is to present a new method 
for the evaluation of the VCA grand potential based on 
the Lanczos method which leads to sufficiently accurate 
results even for larger clusters where full diagonalization 



2 



is no longer possible. Using this method we investigate 
the competing phases in the two-dimensional Hubbard 
model at zero temperature for clusters up to 10 sites. 
This implies a substantial qualitative step forward as 
short-range correlations between different Cooper pairs 
can be included - opposed to calculations based on 2 x 2 
clusters. 



II. VARIATIONAL CLUSTER APPROACH 

The variational cluster approach is one of the possible 
approximation schemes that can be constructed within 
the self-energy-functional theory (SFT).i^ The SFT pro- 
vides a variational scheme to use dynamical information 
from an exactly solvable "reference system" (for exam- 
ple an isolated cluster) to approximate the physics of a 
system in the thermodynamic limit. For a system with 
Hamiltonian H = Ho(t) + Hi(U) and one-particle and 
interaction parameters t and U, the grand potential is 
written as a functional of the self-energy X as 

0[S] = J F[S]+Trln(G(7 1 -S)" 1 , (2) 

with the stationary property (5f2[S p hy S ] = for the phys- 
ical self-energy. Here, Go = (w + \x — t)^ 1 is the free 
Green's function of the original model in the thermo- 
dynamic limit at frequency w, and F[S\ is the Legen- 
dre transform of the universal Luttingcr-Ward functional. 
Due to its universality it is the same as the functional 
for a "simpler" problem with the same interaction but 
a modified one-particle part t' . The stationary solu- 
tions are obtained within the subspace of self-energies 
S = S(t') of the simpler problem that is spanned by 
varying t' . This restriction constitutes the approxima- 
tion. Details of the approach are described in Refs. ll5lfl6l . 

The VCA& is generated within the SFT by choosing as 
a reference system a set of isolated clusters which tile up 
the original infinite lattice. By construction, the VCA 
correctly incorporates correlation effects in the electron 
self-energy up to the length scale given by the cluster size. 
Beyond this scale it acts like a mean-field approximation. 
One of the main advantages of the VCA as compared to 
the simpler cluster perturbation theory^ consists in its 
ability to describe (normal and off-diagonal) long-range 
order by including suitably chosen fictitious symmetry- 
breaking Weiss fields in the set of variational parameters. 
Microscopically coexisting phases can be obtained using 
several Weiss fields. The method links in a consistent way 
the static thermodynamics with the frequency-dependent 
one-particle excitation spectra (photocmission). Details 
of the approach have been described elsewhere.^^ 

The VCA grand potential to be calculated in practice 
reads as 

fi = Q! + Trln (Gq 1 - £) _1 - Trhi(G') . (3) 

Here, Go is the free Green's function of the model given 
by Eq. JT]), fi', S, and G' are the grand potential, the 
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FIG. 1: The SFT grand potential Q of the half-filled (n = 
1) Hubbard model with nearest-neighbor hopping t = — 1 
and U = 8t as a function of the variational parameter /iaf 
(staggered magnetic Weiss field). Reference system (inset): 
L c — 10 clusters. We compare results obtained by integration 
over real frequencies with Lorentzian broadenings r\ = 0.1 
(dotted lines) and r\ = 0.05 (dashed lines), as well as for 
the "Q-matrix" evaluation (see text, solid lines). Sl = 100 
Lanczos iteration steps have been performed. 

sclf-encrgy and the Green's function of the cluster refer- 
ence system which depend on the one-particle parameters 
t' . In the present study we consider clusters with L c = 4, 
8, and 10 sites to search for the stationary solution char- 
acterized by the condition dSl/dt' = O; 10 This stationary 
point provides a good approximation to the exact solu- 
tion for the system in the thcrmodynamical limit if the 
sclf-encrgy is sufficiently "short ranged", i.e. sufficiently 
localized within the cluster. 

As discussed in Refs. Hlsl, it is important to evaluate ft 
with high accuracy in order to resolve the relevant energy 
scales of the competing phases, especially in the electron- 
doped case. Here, we present a method in which this eval- 
uation can be done without a numerical integration over 
frequencies (note that frequency integration is implicit in 
the Trln ■ • ■ terms in Eq. |(3J)). We start from the expres- 
sion derived in Ref . [l5| which expresses the Tr In • • • terms 
as a sum over the single-particle excitation energies: 

Trln (Go 1 - S) _1 = -^Tln(l + e^""*) - R 

T =Y / ^ n e(-UJ m )-B, (4) 
m 

and 

Tr In G' = ln(l + c^"™ ) - R 

m 

T = J2"' m 6(-O-R- (5) 

m 

Here Q(ui) is the Heaviside step function, (3 = 1/T the in- 
verse temperature. uj' m are the one-particle excitation en- 
ergies of the reference system, i.e. the poles of G', and ui m 
are the poles of the VCA Green's function (Gq 1 — S) _1 . 
R represents a contribution due to the poles of the self- 
energy (see Ref. [l5l ) which cancels out in Eq. Q and can 
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thus be ignored. The excitation energies uj' m = E r — E s 
of the reference system (i.e. of the cluster) can be readily 
obtained with the help of the Lanczos algorithm from the 
eigenenergies E r of the reference system. Here, we intro- 
duce the notation m = (r, s), to indicate an excitation 
between two states s and r. The major difficulty consists 
in finding the poles ui m of the VCA Green's function. 

This can be done in the following way: Consider the 
Lchmann representation^ of G which can be cast into 
the formic 



1 



?m/3 



(6) 



where a refers to the one-particle orbitals of the cluster 
(typically a = (site i, spin a) but it can also include an 
orbital index). The "Q-matrix" is defined as: 



T=0 



exp(-/3 J B r ) + exp(-/3£; s ) 



Z' 



<$r,o(0|c Q |s) + S sfi {r\c 



(7) 



The spectral weight (residue) of G' a g(u) at a pole w = u)' m 

is given by Q am Q^ ml3 - Z' = Y. r e ~^ Er is the grand- 
canonical partition function at finite temperature, and |0) 
denotes the (grand-canonical) ground state of the refer- 
ence system. Introducing the diagonal matrix g m n{^) = 
S mn / (w -ijj' m ), we have: 



G\u>) = Qg{uj)Q ] 



(8) 



Defining V = t — t', which in case of the VCA typically 
includes the inter-cluster hopping terms, the "subtrac- 
tion" of the fictitious Weiss fields, as well as shifts of the 
one-particle energies (see below) r- the VCA expression 
for the lattice Green's function can be written as: 



G = 



1 



f 



G„ 1 S 



(G') _1 - V 



(9) 



This expression can be transformed with the help of the 
Q-matrix Eq. © and Eq. ©: 



G = 



1 



(QgQ^) 1 v 

QgQ^ + QgQ) V QgQ} + 

Q{g + g Q ] VQ g+ ) Qt 
1 



Q- 



1 ot 



QWQ 



Note that Q is not a square matrix and that QQ 
-l 



(10) 



Q^Q. Since g 1 = oj — A with A rnn = S mn oj' m , the poles 
of G are now simply given by the eigenvalues of the (fre- 
quency independent) matrix M = A + Q^VQ and can 
be easily found by diagonalization. The dimension of M 
is given by the number of poles of G' with non-vanishing 
spectral weight^ Hence, the above scheme requires the 
knowledge of all excited states of the reference system. 



In Refs. §0, these states have been obtained by a full 
diagonalization of a rather small (2 x 2) cluster. 

For larger clusters, where a full diagonalization is not 
possible, the Lanczos algorithm should, in principle, pro- 
vide precisely the required poles and matrix elements 
Eq. ijTj). In practice, however, there are some difficul- 
ties, as we discuss below. Within the Lanczos method 
the matrix elements G a0 (ui) = ({c a ; Cg)) w of a clus- 
ter Green's function at T = are determined in 2L C 
separate Lanczos procedures^ In each procedure, one 
takes as a Lanczos initial vector one element of the sets 
{ci JO), • • • , c Lc J0)} , {4 JO), • • • , 4 ei JO}} where |0> 

is the cluster ground state. In principle the poles should 
be the same for all matrix elements of the Green's func- 
tion. In practice, however, the poles obtained by the 2L C 
runs are slightly different from each other due to the lim- 
ited numerical accuracy of the Lanczos method. There- 
fore, this kind of Lanczos algorithm is not suited for the 
"Q-matrix" evaluation of the grand potential described 
above, since merging all matrix elements of G' into the 
compact form Eq. ||HJ) would results in a too large matrix 
M that cannot be diagonalizcd. 

Fortunately, the problem can be overcome by means 
of the so-called band Lanczos method^ The difference 
with respect to the standard algorithm is that the sets of 
initial vectors given above are used simultaneously within 
one single Lanczos run. This yields the same set of poles 
for all index pairs (a, j3) as well as the corresponding 
weights. The dimension of the matrix M is given by the 
number of iteration steps 2Sl in the Lanczos procedure. 
In this case, one only needs two Lanczos procedures in- 
stead of 2L C . Using this method, one introduces an error 
due to the limited set (2Sl) of the excited states in the 
reference system that are kept in the Lanczos calculation. 
Generally, however, this error is extremely small since 
excitations with large weight result from states which 
converge very fast with increasing Sl ■ These excitations 
with large weight, on the other hand, are just those which 
are dominant in Eq. ^ compared to excitations with 
small weight.— 

We have checked the accuracy of our method by con- 
sidering the symmetry-broken antiferromagnetic phase 
of the Hubbard model at half-filling (Fig. [T|). One 
can clearly see that the Q-matrix evaluation perfectly 
gives the extrapolation of results obtained by numeri- 
cal frequency integration with finite but small Lorentzian 
broadenings 77. We also verified that the results converge 
very fast with Sl, i-e., typically Sl ~ 100 is fully suf- 
ficient. Last but not least, this improved method sub- 
stantially reduces the computational time. For example, 
a factor of approximately 15 is gained for the L c = 10 
cluster. 



III. RESULTS FOR COMPETING PHASES 

On the basis of this improved evaluation, we investi- 
gate the finite-size behavior of the phase-separation ten- 
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FIG. 4: (Color online) Magnetization M and d-wave order 
parameter D as functions of hole and electron doping for L c — 



FIG. 2: (Color online) Chemical potential fi as function of 
hole doping x. Results for L c = 4 (2 x 2), L c = 8 (4 x 2), 
and L c — 10 clusters. The horizontal dashed lines mark the 
critical /i c , and the vertical dotted lines mark the boundaries 
x\ and X2 of the phase separation region in between. 



TABLE I: Discontinuities Ax, AM, and AD across the PS 
region for hole- and electron doping. 




FIG. 3: (Color online) Same as Fig.[2]but for electron doping. 
Note that there is no phase separation for L c — 10; the dotted 
line marks the quantum critical point. 



dency observed for small clusters in Ref.@. As variational 
parameters we use the Weiss fields /iaf and ft-sc 1° allow 
for antiferromagnetic (AF) and d-wave superconducting 
(SC) orders, respectively^^ as well as an overall shift e of 
the one-particle energies in the cluster to ensure a con- 
sistent treatment of the particle density^ Fig. [5] shows 
our results for the two-dimensional Hubbard model with 
U = 8t and next-nearest-neighbor hopping t' = — 0.3£ for 
the case of hole doping. The calculations have been per- 
formed for L c = 4 (top), L c — 8 (middle), and L c = 10 
clusters (bottom). 

We consider the case of hole-doping first. One can im- 
mediately see that the phase-separation tendency, found 
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for the L c — 4 cluster, weakens progressively when in- 
creasing the cluster size. In particular, the correspond- 
ing "energy scale" A/z = /i* — fj, c diminishes very rapidly 
with increasing L c (A/x = 0.050 for L c = 4, A/j = 0.027 
for L c = 8, and A/x = 0.003 for L c = 10), and appears to 
vanish in the L c — * oo limit. Here, [i* is the point where 
the slope of [x(x) changes sign and fi c the chemical poten- 
tial at the transition point. However, this fact does not 
necessarily imply the absence of macroscopic phase sepa- 
ration in the exact ground state of the model under study. 
As a matter of fact, the exact function /i(x) must have 
a non-positive slope. Therefore, in the phase-separated 
case fJ-(x) becomes a straight line between the two bound- 
aries of the phase-separation region x\ and x^^ Whether 
the exact ground state supports phase separation can be 
derived from the finite-size scaling of the doping disconti- 
nuity Aa; = %?,— x\, see Tab.UJ Unfortunately, no regular 
finite-size behavior can be inferred from Fig.[^]and Tab.U 
for hole doping, probably due to the fact that the clus- 
ters are still too small. Opposed to the clear trend visible 
for the electron-doped case (see Tab. HJ, there is a much 
weaker L c dependence of the discontinuities Aa;, AM, 
and AD, which we rather interpret as being irregular. 
However, our results do not exclude microscopic phase 
separation to persist for L c — > oo. The inclusion of long- 
range Coulomb interaction would then be necessary in 
order to "frustrate" the phase separation occurring in the 
plain Hubbard model and produce microscopic inhomo- 
geneous phases, such as stripes i 23 ! 27 ' 28 We stress that at 
this point only qualitative estimates for L c — > oo rather 
than a convincing finite-size scaling are possible. For a 
discussion on these issues see, e.g., Refs. I24II25II26I . 
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FIG. 5: Evolution of the spectral function upon hole doping 
(L c — 8). Top panel: doping x = 0.025, mixed SC+AF phase. 
Middle panel: x = 0.15, SC phase. Bottom panel: x = 0.25, 
SC phase. A Lorentzian broadening of n = Q.2t has been used 
to display the results. 



The situation is quite different in the electron-doped 
case. Here, not only the phase-separation energy A/z, 
but also the doping discontinuity Ax appears to vanish 
for L c — > 00. In fact, Afj, is already an order of mag- 
nitude smaller than for hole doping in the L c = 4 clus- 
ter.— In addition, already for L c = 10 the transition from 
the AF+SC to the pure SC phase has become continu- 
ous at least within numerical accuracy. In this case, the 
weak phase separation observed at the mean-field level 
for small clusters was simply a signal of a tendency of the 
system to produce microscopically inhomogeneous phases 
(such as stripes), as conjectured in Ref. H. The fact that 
the corresponding energy scale is already very small for 
a small cluster could explain why there is no clear sign 
of stripes in electron-doped materials and could possibly 
be related to the much smaller pseudogap energy scale, 
as discussed in Ref. 

Contrary to the phase-separation energy, the AF and 
SC order parameters M and D plotted in Fig. 2] only dis- 
play a rather weak cluster-size dependence. This shows 
that already a small 2x2 cluster describes the static 
ground-state quantities with a rather good accuracy, ex- 
cept for cases close to a phase transition. Finite-size ef- 
fects are more pronounced for D because the SC order 
parameter is a non-local quantity which converges slower 
with increasing cluster size as compared to M. Never- 
theless, from our results we can argue that for both, hole 
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FIG. 6: Spectral function for x = 0.25 [L c = 8), SC phase, as 
in Fig. [5] but for a different stationary point with a particle 
density n' = 0.75... in the cluster. 



and electron doping, at least substantial SC fluctuations 
remain in the thermodynamic limit even in the AF phase. 
A more precise finite-size scaling to identify as to whether 
one really has long-range SC order is not possible since 
that would involve larger constant cluster shapes (2x2, 
4 x 4,...) which are not accessible by the present VCA. 

The comparison of Fig. Q] with previous calculations^ 
shows that the inclusion of the energy shift e as a vari- 
ational parameter provides results which depend only 
weakly on the cluster size L c . For example, we find a 
mixed AF+SC phase for small doping for all cluster sizes 
considered. This can be understood by the fact that the 
inclusion of additional variational parameters "optimize" 
the ground state of the reference system towards the ex- 
act solution of the infinite lattice. 

Fig. [5] shows the evolution of the single-particle spec- 
tral function upon hole doping calculated with the L c = 8 
reference system. In the upper plot for doping x = 0.025, 
the system is still in the mixed AF+SC phase. This is 
the reason for the "back-turning" of the quasi-particle- 
like band around (tt/2,tt/2) although the chemical po- 
tential already "touches" the band at this wavevector. 
For higher doping (middle and lower panel), we clearly 
see a transition to a dispersion crossing the chemical po- 
tential in the nodal direction, in agreement with angle- 
resolved photoemission experiments. The low-energy co- 
herent quasi-particle band with a width of the order of 
a few times J has been replaced by a band of width of 
a few times t. The qualitative trend is well known from 
QMC calculations^ This represents a clear improvement 
as compared to our previous results for L c = 4 clusters, 
where the dispersion in the nodal direction showed "back- 
turning" signals also for higher dopings. 

In spite of the improvement for the nodal direction, the 
d-wave SC gap appears to be too large for the slightly 
overdoped case x = 0.25. Furthermore, for higher dop- 
ings one expects a decrease of the weight of the upper 
Hubbard band which is much stronger than visible in the 
x = 0.25 spectrum. The reason for these shortcomings is 
probably a too strong admixture of the half-filled cluster 
ground state: In the absence of superconductivity and 
for not too high doping, the particle density of the refer- 
ence system (the isolated cluster) is n' = 1. In our case, 
deviations from cluster half-filling arc introduced due to 
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a non-vanishing SC Weiss field only. For L c = 8 we find 
n! = 0.92 at x = 0.25 which is still close to half-filling. 

A physically better description of the spectral den- 
sity at higher dopings can only be achieved when (in 
the absence of superconductivity) starting from a clus- 
ter ground state with ml < 1. This yields a correspond- 
ing SFT grand potential which has to be compared with 
the SFT grand potential for the n' = 1 stationary point. 
Note that for a vanishing SC Weiss field, the VCA can- 
not give a grand potential that is continuous in the entire 
doping range. This is an artifact of the VCA which levels 
off and eventually becomes irrelevant in the large-cluster 
limit. 

In fact, there is a second stationary point for x = 0.250 
with a particle density in the cluster n' = 0.755. Due to 
the non-vanishing SC Weiss field, this is close but not 
equal to the commensurate cluster filling n 1 = 0.75. This 
stationary point, however, exhibits an SFT grand poten- 
tial which is higher as compared to that of the n' = 0.92 
solution and, consequently, should be disregarded. It is 
nevertheless interesting to discuss the spectral density of 
this (metastable) state which is shown in Fig. [6] As could 
have been expected, the SC gap is much smaller (and ac- 
tually not visible on the scale of the figure) and, as com- 
pared to the corresponding spectrum in Fig. El the weight 
of the upper Hubbard band is clearly reduced. Moreover, 
signatures of magnetic order are no longer visible in this 
spectrum. 



IV. CONCLUSIONS 

We have developed a new method to evaluate the VCA 
grand potential which avoids numerical integrations over 



real or Matsubara frequencies, even for large clusters, for 
which a complete diagonalization is not feasible. This 
provides a sufficient accuracy to study the cluster-size 
dependence of the phase-separation tendency obtained 
in previous works. The results of the present paper sug- 
gest that in the hole-doped case phase separation ob- 
served for small clusters persists for L c — > oo, i.e. in 
the exact ground state, while it eventually disappears in 
case of electron doping. This would explain why there is 
no clear sign of stripes in electron-doped materials, and 
could possibly be related to the much smaller (or even ab- 
sent) pseudogap energy scale with respect to hole-doped 
materials £^ 
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